Opioid-related Deaths 1990-2019

Background on the Opioid Crisis in America

The United States opioid epidemic is a nationwide public health crisis. Initially driven by increased consumption and availability of pharmaceutical opioids such as oxycontin, an increasing number of opioid overdoses are now related to heroin and illicitly manufactured fentanyl and fentanyl analogs. Deaths related to opioid use have sharply increased from 1996 when advertisement around opioid based pharmaceuticals became unrestricted.

Research Questions

I wanted to look at the steady increase in opioid-related deaths in the US, to visualize this I used a simple line graph plotting the data from 1990 to 2019. I also compared opioid related deaths to other drug related deaths in the US along the same time scale, initially just alcohol but then also cocaine and amphetamine. As my final research question I wanted to plot the data spatially by looking at how opioid related deaths have differed across different US states. However, I wanted to keep the temporal aspect of the data, so therefore I made a data file for every year with every states opioid-related death rate

Data Source

All the data was sourced from the Institute for Health Metrics and Evaluation

Data Wrangling

library(here)
library(tidyverse)
library(readr)
library(ggplot2)
library(shiny)

# read and open data
## save data file as df

df <- read.csv(here('rateofdeath.csv'))


#deleting unnessecary columns

columns_to_delete <- c("measure_id", "measure_name","cause_id", "metric_id","cause_name","metric_name","age_id","sex_id","sex_name","age_name","upper","lower")
df <- df[, !(names(df) %in% columns_to_delete)]


#renaming collumns

##location_name change to state
colnames(df)[colnames(df) == "location_name"] <- "state"

##val change to number_deaths 
colnames(df)[colnames(df)=="val"] <-"death_rate"

#change to lower case
df$state <- tolower(df$state)


df <- df %>% mutate(number_deaths = as.numeric(death_rate))


#subset data set by year
##1990
desired_year1990 <- 1990
year1990 <- df[which(df$year == desired_year1990),]

##1991
desired_year1991 <- 1991
year1991 <- df[which(df$year == desired_year1991),]

##1992
desired_year1992 <- 1992
year1992 <- df[which(df$year == desired_year1992),]

##1993
desired_year1993 <- 1993
year1993 <- df[which(df$year == desired_year1993),]

##1994
desired_year1994 <- 1994
year1994 <- df[which(df$year == desired_year1994),]

##1995
desired_year1995 <- 1995
year1995 <- df[which(df$year == desired_year1995),]

##1996
desired_year1996 <- 1996
year1996 <- df[which(df$year == desired_year1996),]

##1997
desired_year1997 <- 1997
year1997 <- df[which(df$year == desired_year1997),]

##1998
desired_year1998 <- 1998
year1998 <- df[which(df$year == desired_year1998),]

##1999
desired_year1999 <- 1999
year1999 <- df[which(df$year == desired_year1999),]

##2000
desired_year2000 <- 2000
year2000 <- df[which(df$year == desired_year2000),]

##2001
desired_year2001 <- 2001
year2001 <- df[which(df$year == desired_year2001),]

##2002
desired_year2002 <- 2002
year2002 <- df[which(df$year == desired_year2002),]

##2003
desired_year2003 <- 2003
year2003 <- df[which(df$year == desired_year2003),]

##2004
desired_year2004 <- 2004
year2004 <- df[which(df$year == desired_year2004),]

##2005
desired_year2005 <- 2005
year2005 <- df[which(df$year == desired_year2005),]

##2006
desired_year2006 <- 2006
year2006 <- df[which(df$year == desired_year2006),]

##2007
desired_year2007 <- 2007
year2007 <- df[which(df$year == desired_year2007),]

##2008
desired_year2008 <- 2008
year2008 <- df[which(df$year == desired_year2008),]

##2009
desired_year2009 <- 2009
year2009 <- df[which(df$year == desired_year2009),]

##2010
desired_year2010 <- 2010
year2010 <- df[which(df$year == desired_year2010),]

##2011
desired_year2011 <- 2011
year2011 <- df[which(df$year == desired_year2011),]

##2012
desired_year2012 <- 2012
year2012 <- df[which(df$year == desired_year2012),]

##2013
desired_year2013 <- 2013
year2013 <- df[which(df$year == desired_year2013),]

##2014
desired_year2014 <- 2014
year2014 <- df[which(df$year == desired_year2014),]

##2015
desired_year2015 <- 2015
year2015 <- df[which(df$year == desired_year2015),]

##2016
desired_year2016 <- 2016
year2016 <- df[which(df$year == desired_year2016),]

##2017
desired_year2017 <- 2017
year2017 <- df[which(df$year == desired_year2017),]

##2018
desired_year2018 <- 2018
year2018 <- df[which(df$year == desired_year2018),]

##2019
desired_year2019 <- 2019
year2019 <- df[which(df$year == desired_year2019),]

Making Maps

# could make maps for every year, however these would all have their individual 
#scales for death rate and therefore the colour changing would only represent change for that specific year and not across the years  

library(maps)
us_states <- map_data("state")

#change region to state in map data
colnames(us_states)[colnames(us_states) == "region"] <- "state"

#create map for 1990
p1990 <- ggplot(data = us_states,
            mapping = aes(x = long, y = lat,
                          group = group))

p1990 + geom_polygon(fill = "white", color = "black")

p1990 <- ggplot(data = us_states,
            aes(x = long, y = lat,
                group = group, fill = state))

p1990 + geom_polygon(color = "gray90", size = 0.1) + guides(fill = FALSE)

#align map to correct latitude and longitude
p1990 <- ggplot(data = us_states,
            mapping = aes(x = long, y = lat,
                          group = group, fill = state))

p1990 + geom_polygon(color = "gray90", size = 0.1) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  guides(fill = FALSE)

#join opioid death data with map data
map1990 <- left_join(year1990, us_states)

#plot opioid death on map data for 1990
p1990 <- ggplot(data = map1990,
            aes(x = long, y = lat,
                group = group, fill = death_rate))

p1990 + geom_polygon(color = "gray90", size = 0.1) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  scale_fill_gradient(low = "white", high = "darkred", na.value = "white", name = "Death Rate per 100,000")

Making maps for each year

#make a loop for all death maps 1990-2019
#keep death rate scale constant through all years to show true change in death rate

overall_scale_limits <- range(df$death_rate, na.rm = TRUE)

generate_and_save_map <- function(current_year) {
  current_data <- df %>% filter(year == current_year)
  map_data <- left_join(current_data, us_states, by = "state")
  
  p <- ggplot(data = map_data,
              aes(x = long, y = lat,
                  group = group, fill = death_rate)) +
    geom_polygon(color = "black", size = 0.1) +
    coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
    scale_fill_gradient(low = "white", high = "darkred", na.value = "yellow", 
                        name = "Death Rate per 100,000", limits = overall_scale_limits) +
    ggtitle(paste("US Opioid-Use Related Death Rate - Year", current_year)) +
    labs(subtitle = "Institute for Health Metrics and Evaluation") +
    theme(plot.title = element_text(size = 18, face = "bold"),
          plot.subtitle = element_text(size = 14),
          panel.grid.major = element_blank(),  # Remove major gridlines
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "transparent"))
  
  # Save the plot as an image
  ggsave(filename = paste("opioid_death_map_", current_year, ".png", sep = ""), plot = p)
}

# Generate and save frames for each year
lapply(1990:2019, generate_and_save_map)

## [[1]]
## [1] "opioid_death_map_1990.png"
## 
## [[2]]
## [1] "opioid_death_map_1991.png"
## 
## [[3]]
## [1] "opioid_death_map_1992.png"
## 
## [[4]]
## [1] "opioid_death_map_1993.png"
## 
## [[5]]
## [1] "opioid_death_map_1994.png"
## 
## [[6]]
## [1] "opioid_death_map_1995.png"
## 
## [[7]]
## [1] "opioid_death_map_1996.png"
## 
## [[8]]
## [1] "opioid_death_map_1997.png"
## 
## [[9]]
## [1] "opioid_death_map_1998.png"
## 
## [[10]]
## [1] "opioid_death_map_1999.png"
## 
## [[11]]
## [1] "opioid_death_map_2000.png"
## 
## [[12]]
## [1] "opioid_death_map_2001.png"
## 
## [[13]]
## [1] "opioid_death_map_2002.png"
## 
## [[14]]
## [1] "opioid_death_map_2003.png"
## 
## [[15]]
## [1] "opioid_death_map_2004.png"
## 
## [[16]]
## [1] "opioid_death_map_2005.png"
## 
## [[17]]
## [1] "opioid_death_map_2006.png"
## 
## [[18]]
## [1] "opioid_death_map_2007.png"
## 
## [[19]]
## [1] "opioid_death_map_2008.png"
## 
## [[20]]
## [1] "opioid_death_map_2009.png"
## 
## [[21]]
## [1] "opioid_death_map_2010.png"
## 
## [[22]]
## [1] "opioid_death_map_2011.png"
## 
## [[23]]
## [1] "opioid_death_map_2012.png"
## 
## [[24]]
## [1] "opioid_death_map_2013.png"
## 
## [[25]]
## [1] "opioid_death_map_2014.png"
## 
## [[26]]
## [1] "opioid_death_map_2015.png"
## 
## [[27]]
## [1] "opioid_death_map_2016.png"
## 
## [[28]]
## [1] "opioid_death_map_2017.png"
## 
## [[29]]
## [1] "opioid_death_map_2018.png"
## 
## [[30]]
## [1] "opioid_death_map_2019.png"

#Making a gif

# Install and load the magick package
library(magick)

# Get the list of PNG files in the current directory
file_names <- list.files(pattern = "\\.png$", full.names = TRUE)

# Read the PNG files into R as magick images
images <- image_read(file_names)


animation <- image_animate(images, fps = 2)

# Specify the output format as GIF when writing the animation
animation_file <- "output.gif"
image_write(animation, animation_file, format = "gif")

Making an Interactive Map

#making maps that change by year and you can look at inidividual death rate for each state

library(shiny)
library(shinyWidgets)

map_changing <- fluidPage(
  titlePanel("US Opioid-Use Related Death Rate"),
  
  # Create a slider input for selecting the year
  sliderTextInput("year", "Select Year:",
                  choices = as.character(1990:2019), selected = "1990"),
  
  # Display the plot
  plotOutput("map")
)

server <- function(input, output) {
  # Generate the plot based on the selected year
  output$map <- renderPlot({
    current_data <- df %>% filter(year == input$year)
    map_data <- left_join(current_data, us_states, by = "state")
    
    ggplot(data = map_data,
           aes(x = long, y = lat,
               group = group, fill = death_rate)) +
      geom_polygon(color = "black", size = 0.1) +
      coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
      scale_fill_gradient(low = "white", high = "red", na.value = "yellow", 
                          name = "Death Rate", limits = overall_scale_limits) +
      ggtitle(paste("US Opioid-Use Related Death Rate - Year", input$year)) +
      labs(subtitle = "Institute for Health Metrics and Evaluation") +
      theme(plot.title = element_text(size = 18, face = "bold"),
            plot.subtitle = element_text(size = 14),
            panel.grid.major = element_blank(),  
            panel.grid.minor = element_blank(),
            panel.background = element_rect(fill = "transparent"))
  })
}

shinyApp(map_changing, server)

Shiny applications not supported in static R Markdown documents

#Combine shiny and plotly to make an interactive map where you can hover over the states to see a value and also go through year by year



library(plotly)

map_changing <- fluidPage(
  titlePanel("US Opioid-Use Related Death Rate"),
  
  # Create a slider input for selecting the year
  sliderTextInput("year", "Select Year:",
                  choices = as.character(1990:2019), selected = "1990"),
  
  # Display the plot
  plotlyOutput("map")
)

server <- function(input, output) {
  # Generate the plot based on the selected year
  output$map <- renderPlotly({
    current_data <- df %>% filter(year == input$year)
    map_data <- left_join(current_data, us_states, by = "state")
    
    p <- ggplot(data = map_data,
                aes(x = long, y = lat,
                    group = group, fill = death_rate, text = paste("State: ", state, "<br>Death Rate: ", death_rate))) +
      geom_polygon(color = "black", size = 0.1) +
      coord_map() +
      scale_fill_gradient(low = "white", high = "darkred", na.value = "yellow", 
                          name = "Death Rate", limits = overall_scale_limits) +
      ggtitle(paste("US Opioid-Use Related Death Rate - Year", input$year)) +
      labs(subtitle = "Institute for Health Metrics and Evaluation") +
      theme(plot.title = element_text(size = 18, face = "bold"),
            plot.subtitle = element_text(size = 14),
            panel.grid.major = element_blank(),  
            panel.grid.minor = element_blank(),
            panel.background = element_rect(fill = "transparent"))
    
    ggplotly(p) %>% layout(geo = list(projection = list(type = 'albers usa'))) 
  })
}



shinyApp(map_changing, server)

Shiny applications not supported in static R Markdown documents